!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Interface to the quadrature formulas.
!!
!! \n
!!
!! This module provides an interface to the collection of quadrature
!! formulas for the computation of \f$d\f$-dimensional integrals. Due
!! to the difficulty of the computations of accurate quadrature
!! formulas, this module relies on a file based database, except for
!! the one-dimensional formulas, which computation "on the fly" is
!! straightforward. All the considered quadrature formulas are
!! symmetric, which makes them suitable for element integrals as well
!! as for side integrals, and the database also includes permutation
!! tables. A permutation table is an integer array which establish a
!! correspondence between vertex permutations, taken in lexicographic
!! order, and quadrature nodes permutations (see the source code for
!! an example). All the quadrature formulas are referred to the
!! canonical simplex (see mod_master_el).
!<----------------------------------------------------------------------
!
! To clarify the concept of a permutation table, let us consider te
! following example: a 4 point quadrature rule on a triangle, with the
! following ordering of vertices (numbers) and quadrature nodes
! (letters):
!
!              3
!              /\
!             /D \
!            /    \
!           /      \
!          /        \
!         /    A     \
!        /            \
!       /_B__________C_\
!      1                2
!
! The vertex permutation to node permutation table is
!      1,2,3  ->  A,B,C,D
!      1,3,2  ->  A,B,D,C
!      2,1,3  ->  A,C,B,D
!      2,3,1  ->  A,C,D,B
!      3,1,2  ->  A,D,B,C
!      3,2,1  ->  A,D,C,B
!-----------------------------------------------------------------------
module mod_numquad

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   read_octave_al

 use mod_linal, only: &
   mod_linal_initialized, &
   eig,   &
   sort

 use mod_perms, only: &
   mod_perms_initialized, &
   fact

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_numquad_constructor, &
   mod_numquad_destructor,  &
   mod_numquad_initialized, &
   get_quad_nodes, &
   get1d_gauss_nodes, get1d_gausslobatto_nodes

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! private members
 !> This is the root directory of the database of the quadrature
 !! formulas. The full structure of the database is
 !! <tt>database_dir/<nd>D/Q<k>.octxt</tt>, where \c nd is the space
 !! dimension: 2, 3, \f$\ldots\f$ and \c k is the exactness of the
 !! rule (with two digits).
 !!
 !! Each quadrature formula contains three variables:
 !! <ul>
 !!  <li>\c x: the quadrature nodes in a <tt>(d,m)</tt> real array,
 !!  where \c d is the spatial dimension and \c m is the number of
 !!  nodes
 !!  <li>\c w: the quadrature weights in a <tt>(m)</tt> real array
 !!  <li>\c stab: the permutation table in a <tt>(\f$(d+1)!\f$,m)</tt>
 !!  integer array
 !! </ul>
 !<
 character(len=*), parameter :: &
    database_dir = './QUADRATURE'

! Module variables

 ! public members
 logical, protected ::               &
   mod_numquad_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_numquad'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_numquad_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)  .or. &
       (mod_kinds_initialized.eqv..false.)     .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_octave_io_initialized.eqv..false.) .or. &
       (mod_linal_initialized.eqv..false.)     .or. &
       (mod_perms_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_numquad_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_numquad_initialized = .true.
 end subroutine mod_numquad_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_numquad_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_numquad_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_numquad_initialized = .false.
 end subroutine mod_numquad_destructor

!-----------------------------------------------------------------------
 
 !> Return a quadrature rule for integration on the
 !! \f$d\f$-dimensional canonical simplex with exactness deg.
 !<
 subroutine get_quad_nodes(d,deg,x,w,stab)
  integer, intent(in) :: d, deg
  real(wp), intent(out), allocatable :: x(:,:), w(:)
  integer, intent(out), allocatable, optional :: stab(:,:)
 
  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'get_quad_nodes'
  character(len=len(database_dir)+1+1+3+2+6) :: file_name
  character(len=100) :: message

   ! The cases d=0 and d=1 are simple; for d.ge.2, the lowest order
   ! cases can still be dealt with on the fly; everything else must be
   ! read from a database of precomputed quadrature formulas.
   if(d.eq.0) then

     ! We define one quadrature point with no coordinates. This has
     ! two advantages:
     ! - it is consistent with the definition of the affine mapping
     !   given in mod_grid, which then reduces to a constant
     ! - it indicates that there is one quadrature point, thus
     !   permitting the evaluation of the boudary integral as a
     !   summation on the quadrature points without any special
     !   casing.
     allocate( x(0,1) , w(1) ); w = 1.0_wp
     if(present(stab)) then
       allocate( stab(1,1) ); stab = 1
     endif

   elseif(d.eq.1) then

     call get1d_gauss_nodes(deg,x,w,stab)

   elseif(deg.le.1) then ! barycenter

     allocate( x(d,1) , w(1) )
     x = 1.0_wp/real(d+1,wp)
     w = 1.0_wp/real(fact(d),wp)
     if(present(stab)) then
       allocate( stab(fact(d+1),1) ); stab = 1
     endif

   elseif((d.eq.2).and.(deg.le.3)) then

     d2_special_cases: select case(deg)
      case(2)
       allocate( x(2,3) , w(3) )
       x(1,:) = (/ 1.0_wp/6.0_wp , 2.0_wp/3.0_wp , 1.0_wp/6.0_wp /)
       x(2,:) = (/ 1.0_wp/6.0_wp , 1.0_wp/6.0_wp , 2.0_wp/3.0_wp /)
       w = 0.5_wp/3.0_wp
       if(present(stab)) then
         allocate( stab(6,3) )
         stab(1,:) = (/ 1 , 2 , 3 /)
         stab(2,:) = (/ 1 , 3 , 2 /)
         stab(3,:) = (/ 2 , 1 , 3 /)
         stab(4,:) = (/ 2 , 3 , 1 /)
         stab(5,:) = (/ 3 , 1 , 2 /)
         stab(6,:) = (/ 3 , 2 , 1 /)
       endif
      case(3)
       allocate( x(2,4) , w(4) )
       x(1,:) = (/ 1.0_wp/3.0_wp , &
                   1.0_wp/5.0_wp , 3.0_wp/5.0_wp , 1.0_wp/5.0_wp /)
       x(2,:) = (/ 1.0_wp/3.0_wp , &
                   1.0_wp/5.0_wp , 1.0_wp/5.0_wp , 3.0_wp/5.0_wp /)
       w = 0.5_wp*(/ -9.0_wp/16.0_wp , &
                     25.0_wp/48.0_wp*(/ 1.0_wp , 1.0_wp , 1.0_wp /) /)
       if(present(stab)) then
         allocate( stab(6,4) )
         stab(1,:) = (/ 1 , 2 , 3 , 4 /)
         stab(2,:) = (/ 1 , 2 , 4 , 3 /)
         stab(3,:) = (/ 1 , 3 , 2 , 4 /)
         stab(4,:) = (/ 1 , 3 , 4 , 2 /)
         stab(5,:) = (/ 1 , 4 , 2 , 3 /)
         stab(6,:) = (/ 1 , 4 , 3 , 2 /)
       endif
     end select d2_special_cases

   else

     write(file_name,'(a,a,i1,a,i2.2,a)') &
               database_dir,'/',d,'D/Q',deg,'.octxt'
     call new_file_unit(fu,ierr)
     open(fu,file=trim(file_name), &
          status='old',action='read',form='formatted',iostat=ierr)
     if(ierr.ne.0) then
       write(message,'(A,A,A)') &
         'Problem opening file "',trim(file_name),'".'
       call error(this_sub_name,this_mod_name,message)
     endif
     call read_octave_al(x,'x',fu)
     call read_octave_al(w,'w',fu)
     if(present(stab)) call read_octave_al(stab,'stab',fu)
     close(fu,iostat=ierr)
     if(ierr.ne.0) then
       write(message,'(A,A,A)') &
         'Problem closing file "',trim(file_name),'".'
       call error(this_sub_name,this_mod_name,message)
     endif

   endif

   if(minval(w).lt.0.0_wp) then
     write(message,'(A,I3,A)') 'Quadrature formula for deg=',deg, &
       ' contains negative weights.'
     call warning(this_sub_name,this_mod_name,message)
   endif

 end subroutine get_quad_nodes

!-----------------------------------------------------------------------
 
 !> Gauss quadrature nodes.
 !!
 !! See <a href="http://dx.doi.org/10.1137/1015032">[Golub, 1973]</a>
 !! and <a href="http://dx.doi.org/10.1145/174603.174605">[Gautschi,
 !! 1994]</a> for detailed explanations of the algorithm.
 subroutine get1d_gauss_nodes(deg,x,w,stab)
  integer, intent(in) :: deg
  real(wp), intent(out), allocatable :: x(:,:), w(:)
  integer, intent(out), allocatable, optional :: stab(:,:)
 
  integer :: n, ierr, i
  integer, allocatable :: ii(:)
  real(wp), allocatable :: a(:), b(:), jacm(:,:), wi(:), &
    zr(:,:), zi(:,:)

   n = ceiling(0.5_wp*real(deg+1,wp)) ! number of nodes
   allocate(x(1,n),w(n))
   ! The nodes are first compute on [-1,1]
   if(n.eq.1) then
     x = 0.0_wp
     w = 2.0_wp
   else
     call coeflege(n,a,b)
     ! notice that x and w are already allocated
     allocate(jacm(n,n),wi(n),zr(n,n),zi(n,n))
     allocate(ii(n))
     jacm = 0.0_wp
     ! set the three diagonals
     jacm(1,1) = a(1)
     jacm(1,2) = sqrt(b(2))
     do i=2,n-1
       jacm(i,i-1) = sqrt(b(i))
       jacm(i,i)   = a(i)
       jacm(i,i+1) = sqrt(b(i+1))
     enddo
     jacm(n,n-1) = sqrt(b(n))
     jacm(n,n) = a(n)
     call eig(jacm,x(1,:),wi,zr,zi,ierr)
     ii = (/ (i, i=1,n) /)
     call sort(x(1,:),ii)
     w = 2.0_wp*(zr(1,ii)**2)
   endif

   ! Shift the interval
   x = 0.5_wp*x + 0.5_wp
   w = 0.5_wp*w
 
   if(present(stab)) then
     allocate(stab(2,n))
     stab(1,:) = (/ (i, i=1,n,+1) /)
     stab(2,:) = (/ (i, i=n,1,-1) /)
   endif
 end subroutine get1d_gauss_nodes
 
!-----------------------------------------------------------------------
 
 pure subroutine coeflege(n,a,b)
 ! coefficients of the Legendre polynomial
  integer, intent(in) :: n
  real(wp), intent(out), allocatable :: a(:), b(:)

  integer :: k
 
  allocate(a(n),b(n))
  a = 0.0_wp
  b(1) = 2.0_wp
  do k=2,n
    b(k) = 1.0_wp/(4.0_wp - 1.0_wp/(real(k-1,wp)**2))
  enddo

 end subroutine coeflege

!-----------------------------------------------------------------------
 
 !> Gauss-Lobatto nodes
 !!
 !! \c deg is the degree of exactness of the formula, which is equal
 !! to \f$2n-3\f$, where \f$n\f$ is the number of points. In practice,
 !! the number of points is computed as
 !! \f{displaymath}{
 !!  n = {\rm ceiling}\left(\frac{{\rm deg}+3}{2}\right).
 !! \f}
 pure subroutine get1d_gausslobatto_nodes(deg,x,w,stab)
  integer, intent(in) :: deg
  real(wp), intent(out), allocatable :: x(:,:), w(:)
  integer, intent(out), allocatable, optional :: stab(:,:)
 
  integer :: n, ierr, i
  integer, allocatable :: ii(:)
  real(wp) :: p0l, p0r, p1l, p1r, pm1l, pm1r, det
  real(wp), allocatable :: a(:), b(:), jacm(:,:), wi(:), &
    zr(:,:), zi(:,:)
  character(len=*), parameter :: &
    this_sub_name = 'get1d_gausslobatto_nodes'
 
   n = ceiling(0.5_wp*real(deg+3,wp)) ! number of nodes
   allocate(x(1,n),w(n))
   ! The nodes are first compute on [-1,1]
   if(n.eq.1) then
     x = 0.0_wp ! this is not really Gauss-Lobatto, but anyway...
     w = 2.0_wp
   elseif(n.eq.2) then
     x(1,:) = (/ -1.0_wp , 1.0_wp /)
     w      = (/  1.0_wp , 1.0_wp /)
   else
     call coeflege(n-1,a,b) ! first coefficients as in Gauss
     ! Here we need to compute two additional coefficients
     p0l = 0.0_wp; p0r = 0.0_wp; p1l = 1.0_wp; p1r = 1.0_wp
     do i=1,n-1
       pm1l = p0l; p0l  = p1l
       pm1r = p0r; p0r  = p1r
       p1l = (-1.0_wp-a(i))*p0l - b(i)*pm1l
       p1r = ( 1.0_wp-a(i))*p0r - b(i)*pm1r
     enddo
     det = p1l*p0r-p1r*p0l
     ! wi is used as temporary to reallocate a,b
     allocate(wi(size(a))); wi = a
     deallocate(a); allocate(a(size(wi)+1))
     a = (/ wi , (-1.0_wp*p1l*p0r-1.0_wp*p1r*p0l)/det /)
     wi = b; deallocate(b); allocate(b(size(wi)+1))
     b = (/ wi ,         2.0_wp*p1l*p1r/det           /)
     deallocate(wi)

     ! Now everything proceeds as for the Gauss nodes
     allocate(jacm(n,n),wi(n),zr(n,n),zi(n,n))
     allocate(ii(n))
     jacm = 0.0_wp
     ! set the three diagonals
     jacm(1,1) = a(1)
     jacm(1,2) = sqrt(b(2))
     do i=2,n-1
       jacm(i,i-1) = sqrt(b(i))
       jacm(i,i)   = a(i)
       jacm(i,i+1) = sqrt(b(i+1))
     enddo
     jacm(n,n-1) = sqrt(b(n))
     jacm(n,n) = a(n)
     call eig(jacm,x(1,:),wi,zr,zi,ierr)
     ii = (/ (i, i=1,n) /)
     call sort(x(1,:),ii)
     w = 2.0_wp*(zr(1,ii)**2)
   endif

   ! Shift the interval
   x = 0.5_wp*x + 0.5_wp
   w = 0.5_wp*w
 
   if(present(stab)) then
     allocate(stab(2,n))
     stab(1,:) = (/ (i, i=1,n,+1) /)
     stab(2,:) = (/ (i, i=n,1,-1) /)
   endif

 end subroutine get1d_gausslobatto_nodes
 
!-----------------------------------------------------------------------

end module mod_numquad

